Site Dilution in the Half-Filled One-Band Hubbard Model: Ring Exchange, Charge 
Fluctuations and Application to La 2 Cui_ a; (Mg/Zn) a .04 
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We study the ground state quantum spin fluctuations around the Neel ordered state for the one- 
band (t,U) Hubbard model on a site-diluted square lattice. An effective spin Hamiltonian, Hs 4 \ is 
generated using the canonical transformation method, expanding to order t(t/U) 3 . Hs^ contains 
four-spin ring exchange terms as well as second and third neighbor bilinear spin-spin interactions. 
Transverse spin fluctuations are calculated to order 1/5 using a numerical real space algorithm 
first introduced by Walker and Walsteadt. Additional quantum charge fluctuations appear to this 
order in t/U, coming from electronic hopping and virtual excitations to doubly occupied sites. 
The ground state staggered magnetization on the percolating cluster decreases with site dilution x, 
vanishing very close to the percolation threshold. We compare our results in the Heisenberg limit, 
t/U — > 0, with quantum Monte Carlo (QMC) results on the same model and confirm the existence of 
a systematic rr-dependent difference between 1/5 and QMC results away from x — 0. For finite t/U, 
we show that the effects of both the ring exchange and charge fluctuations die away rapidly with 
increasing t/U. We use our finite t/U results to make a comparison with results from experiments 
on La 2 Cui_ :c (Mg/Zn) :c 04. 



I. INTRODUCTION 

A. Random disordered magnets 

Magnetic materials and model magnetic systems are 
perhaps the best test benches for the study of collec- 
tive phenomena in nature. This is particularly true in 
the context of systems with frozen or quenched random 
disordepi^. Here, questions such as the sharpness of 
phase transitions in disordered systems 3 -, the stability 
of ground-state symmetry-breaking (random-field) per- 
turbations^ and spin glass behavior arising from random 
frustration^ have come under sharp scrutiny over the past 
thirty years. 

The 1987 discovery of high-temperature superconduc- 
tivity in doped antifcrromagnctic copper oxide materials 
generated a huge amount of interest in quantum anti- 
ferromagnets which continues to this day££ Here, the 
magnetic properties depend strongly on the different pos- 
sible types of quenched disorder and this has proven to 
be a rich source of novel quantum phenomena. An im- 
portant area of investigation has been to explore how the 
ground state of insulating quantum magnets evolves as 
the level of random disorder is changed. The following 
examples represent a small subset of this class of studies. 
A large effort has been targeted towards understanding 
the properties of antiferromagnetic spins chains subject 
to various types of disorde r 8 ^ 10 . Further work investi- 
gated how long range order develops in two and three 
dimensional arrays of weakly coupled integer (Haldanc) 
spin chains^ and even-leg ladders 1 - 1 -. The question of how 
Nccl order develops upon magnetically diluting pure sys- 



tems with quantum spin liquid ground states is another 
field of intensive studyi 2 -. 

Theoretical problems relating to various types of ran- 
dom bond disorder, as opposed to the more material- 
relevant case of site dilution, have also been investi- 
gate d 13 ' 14 . In three dimensional systems, one note- 
worthy example is the so-called antiglass phenomenon 
in LiHo a; Yi_ 2 ;F4 where, for a low concentration, x, of 
magnetic Ho 3+ ions, the dipolar spin glass phase seem- 
ingly disappears 1 ^. Another interesting problem con- 
cerns the role frozen random impurities may play at con- 
ventional and deconfined quantum critical points in two- 
dimensional antifcrromagncts 1 ^. 

However, among the multitude of interesting problems, 
a particular one, possibly because of its seemingly sim- 
ple physical setting and its broad conceptual appeal, has 
drawn considerable attention: that of the evolution of the 
antiferromagnetic Neel ground state in the site-diluted 
S = 1/2 nearest neighbor square lattice Heisenberg anti- 
ferromagnet (SLHAF). 



B. Site-diluted SLHAF and La 2 Cui ^ (Mg/Zn)^ 

As the insulating and antifcrromagnctic parent of 
high-temperature superconductivity in La2- a: Sr a ;Cu04, 
La2Cu04 has quasi two dimensional magnetic exchange 
interactions and a good starting point for its descrip- 
tion is to treat the CuC>2 planes as decoupled SLHAFs. 
Hence, experimental studies on zinc and magnesium sub- 
stitution for copper in I^CuO V > 18 , provided some of 
the earliest motivation and interest in the problem of 
site-diluted SLHAFs 1 -. In particular, Chcong et alr^ 
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found from bulk thermodynamic measurements that in 
the diluted 5=1/2 quantum antiferromagnetic mate- 
rials. La2Cui_ T Zn :E 04 and La2Cui_ 2; Mg a ;04, the Neel 
temperature, Tn, vanishes faster than in other materi- 
als that can be considered as site-diluted classical square 
lattice magnetic systems (either because they have large 
spin 5, or because they have large Ising anisotropics). 

Most importantly, these early experimental results sug- 
gested that Tn, hence long range antiferromagnetic Neel 
order, may vanish at a critical impurity concentration x c 
less than the site dilution percolation threshold for the 
square lattice, x p « 0.41. This possibility was seemingly 
supported by subsequent muon spin relaxation (/iSR) 
and nuclear quadrupole resonance (NQR) experiments 2 ^, 
with these latter measurements also suggesting the possi- 
bility of a second transition below T^(x) into a spin-glass 
like state. 

From a classical point of view, the ground state of 
the SLHAF has two-sublattice Neel order for all x < 
Consequently, early experiments on site-diluted 



La2Cu04 from Refs. 18lf20| implied that either a novel 
quantum ground state develops in the site-diluted SL- 
HAF for x c < x < x p , or that frustrating further neigh- 
bor exchange interactions are important in the real ma- 
terial and that these drive the system into a two di- 
mensional Hcisenbcrg spin glass ground state, presum- 
ably via the proliferation of Villain canted states for 
x c < x < Xj ^^ A ^ e 

The idea that Neel order could disappear in the di- 
luted SLHAF, due to quantum effects, for a concentra- 
tion of magnetic moments less than the geometric site 
percolation threshold (x < x p ) had been suggested by 
some^l, but not all^, early calculations. However, in 
strong contrast to the early body of experimental evi- 
denc o 18 ' 20 and theoretical suggestions^ 7 -, large scale quan- 
tum Monte Carlo (QMC) simulations on the diluted SL- 
HAF find that Neel order survives up to the percolation 
threshold X p 28 ^ 29 . Further, contrary to earlier experi- 
ment a 18 ' 20 , recent neutron scattering studies on single 
crystals of La2Cui_ 2: (Mg/Zn) :c 04 found that long-range 
Neel order does survive up to at least x = 0.39, if not up 
to x p 30 ' 31 . Interestingly, recent QMC studies show that 
the same scenario holds for homogeneous bond dilution, 
with exotic quantum phases appearing only for inhomo- 
gencous dilution where local ladder structures formi^. 

A proposed explanation for the discrepancy between 
the earlier experimcnt o 17 ' 18 ' 20 and the more recent 
one o 30 i 31 is that samples are extremely sensitive to excess 
oxygen, or off-stoichiomctric S, La2Cu04+5, as Cu 2+ is 
substituted by either Zn 2+ or Mg 2+ . Off-stoichiometry 
with S > is hole-doping, which is extremely detrimental 
to long range Neel order. Thus, the present picture, sup- 
ported by both numerica l 19 ' 28 ! 29 and experimenta l 30 ' 31 
studies, is that Neel order survives in the site-diluted SL- 
HAF 2 ^ and La 2 Cui_ :c (Mg/Zn) :r O 4 20»2I up to x p , with 
no intervening exotic quantum ground state for x < x p . 



C. Quantum Monte Carlo vs La2Cui (Mg/Zn^O^ 



While both high precision quantum Monte Carlo 
(QMC) studies of the site-diluted SLHAE22 and neutron 
scattering experiment o 30 ' 31 on LaaCui-^Mg/Zn^C^ 
now find that the Neel order survives up to x c ~ x p (ex- 
actly x c = x p for the QMC simulations), the quantitative 
agreement stops here. There is a systematic discrepancy 
between QMC and the neutron results for the sublattice 
Neel order parameter, [M(x)], as a function of x. The ex- 
perimental and numerical data are reproduced in Fig. [T] 
In this figure, the QMC results of Ref. [2!| are shown by 
the upper solid line. The experimental results (neutron, 
squares, from Ref. [3(j; NQR, triangle, from Ref. [2(J) lie 
on the dashed line, which is a guide to the eye parameter- 
ized by [M(x)]/M(0) = (1 - x/x P Y" a . The QMC results 
lie above the experimental data over the whole range 
< x < x p , as illustrated by the shaded region. Tak- 
ing it as a premise that the QMC data are essentially the 
exact results for the diluted 5 = 1/2 SLHAF, the system- 
atic difference between them and the experimental data 
shown in Fig. [1] suggests that Zn 2+ and Mg 2+ substituted 
Cu 2+ in La2Cu04 are not quantitatively described by 
a site-diluted nearest neighbor Hcisenbcrg Hamiltonian. 
The nature of the discrepancy is in itself interesting. It 
is initially small at low x, increases and reaches a max- 
imum for x ~ 0.35, and decreases upon approaching x p 
such that the "true" underlying microscopic Hamiltonian 
describing La2Cui_ a; (Mg/Zn) ;E 04 seem to also possess a 
percolation threshold very close to that of the idealized 
SLHAF. 
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FIG. 1: Ground state staggered magnetization, [M(x)\, 
as a function of concentration, x, of Zn and Mg in 
La2Cui_ I (Mg/Zn) I 04, normalized to the value for zero dilu- 
tion, M(0)— . The solid line shows the results from quantum 
Monte Cartel for the site-diluted SLHAF. The fi gure is re- 
produced from Ref. [30( ] . 
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D. Ring exchange interactions 

One class of candidate perturbations that may be giv- 
ing the missing physics of diluted La2Cu04 are ring, 
or cyclic, exchange interactions involving multiple inter- 
actions around closed plaquettes of the square lattice. 
Such interactions have received intensive attention re- 
centl y 32 ' 33 ' 34 and have been shown to play an impor- 
tant role in the quantitative description of undiluted 
La2Cu0 4 35 i 36 ' 37 i 38 . Taking as a starting point the one- 
band half-filled Hubbard mode l 39 ' 40 ' 41 , the lowest order 
ring exchange interaction takes its origin in virtual elec- 
tronic hopping process, fourth order in t/U, taking elec- 
trons coherently around a closed square plaquette. Here 
t is the nearest neighbor hopping constant and U is the 
on-site Coulomb energy. Taking it as plausibl o 36 ' 37 ' 38 
that ring exchange is indeed present and a leading per- 
turbation beyond the Hciscnberg model description of 
La2Cu04, it is natural to ask what its effect is on the 
Neel order parameter upon substituting Cu by a concen- 
tration x of non-magnetic ions (see Fig. [1} . This ques- 
tion, which to the best of our knowledge has so far not 
been investigated, is the one that we explore in this pa- 
per. 

To tackle this question, one must return to a problem 
of correlated electrons. The reason is that the spin-only 
Hamiltonian with ring exchange derives from a set of 
electronic hops. As we show below, the elimination of an 
intermediate site in an electron hopping pathway, affects 
the resulting effective spin Hamiltonian in a nontrivial 
manner. Specifically, we consider the problem of a site- 
diluted half-filled one-band Hubbard model away from 
the Heisenberg t/U — * limit. Since here ring exchange 
originates solely from correlated nearest neighbor elec- 
tronic hops, they cannot move the percolation threshold 
to a larger value than the nearest neighbor threshold x p . 
From this constraint alone, ring exchange is an admissi- 
ble candidate for a perturbation to the diluted 5=1/2 
SLHAF, as it preserves the same geometric percolation 
threshold the nearest neighbor Heisenberg model. 

The presence of the ring exchange and second and third 
nearest neighbor bilinear exchange terms in the Hamilto- 
nian, generated by hopping processes to fourth order in 
t/U, leads to a sign problem for currently available QMC 
methods using the standard S z basis representation of 
the Hamiltonian^ 2 -. A direct attack on the site-diluted 
ring exchange Hamiltonian via QMC, such as done for 
the site-diluted Heisenberg model^ is therefore not pos- 
sible at this time. As a first step in investigating the 
role played by ring exchange in the site-diluted Hubbard 
model, we carry out a finite-lattice spin wave calculation 
to order 1/5 on an extended effective spin Hamiltonian 
generated from up to four hop electronic pathways. To 
proceed, we use a real space linear spin wave method 
adapted to finite-size diluted lattices, first developed by 
Walker and Walsteadt^ 3 - in the context of spin glasses and 
similar to that used for the site-diluted nearest neighbor 
Heisenberg antifcrromagnct on the squared and honey- 



comb lattices^. We investigate the role of ring exchange 
on the dependence of the ground state staggered mag- 
netization, [M(a;)], as a function of x. In Ref. [44j], it 
was found that there is a systematic difference between 
the value of this quantity, calculated via the spin wave 
method and the essentially exact QMC- 9 -. From this, it is 
clear that a similar systematic difference should also ex- 
ist between our data calculated using a 1/5 expansion, 
and what would be the not yet available numerically ex- 
act value for [M(x)], as a function of dilution, for the 
extended Hamiltonian. Hence, although the main mo- 
tivation for this project comes from the experiment on 
La2Cui_ 2 -(Mg/Zn) a; 04 — , some care has to be taken in 
attempting to make a direct comparison with experimen- 
tal results. Rather, our results for the extended Hamilto- 
nian and electronic hopping, can be quantitatively bench- 
marked by a comparison with those for the site-diluted 
Heisenberg model, using the same real space expansion 
technique. From a broader perspective, our work pro- 
vides a first glimpse at the role of charge correlation ef- 
fects in the problem of diamagnetic site dilution in the 
one-band Hubbard model. 



E. Charge fluctuations 

The generation of ring and further neighbor exchange 
interactions is not the only effect of extending the anal- 
ysis of the one-band Hubbard model beyond the Heisen- 
berg limit using a perturbation expansion in t/U. We 
have previously shown that extending the expansion to 
order (t/U) 4 generates quantum charge fluctuations 4 ^ 
that are independent of the transverse spin fluctuations of 
localized 5=1/2 moments. These fluctuations appear in 
the perturbation expansion on the square lattice because, 
to this order, the ground state wave function contains 
an admixing with excited states corresponding to doubly 
occupied sites. As doubly occupied sites carry no mo- 
ment, the expectation value for the magnetic moment of 
the Hubbard model is reduced below that expected from 
the effective spin-only Hamiltonian describing transverse 
spin fluctuations. We show here that these charge fluctu- 
ations are a key element in the ultimate success of com- 
parisons between the one-band Hubbard model and ex- 
periments on both undiluted and site-diluted La2Cu04. 
Just as for ring exchange effects, we find that the effects 
of charge fluctuations disappear as the site percolation 
threshold is approached, as four hop electronic processes 
are interrupted by the dilution well before this limit is 
reached. 

The rest of the paper is organized as follows: before 
launching into the calculations, we discuss in Sections 
IIA and IIB some of the caveats that arise when consider- 
ing a low energy effective spin-only Hamiltonian derived 
from a site-diluted Hubbard model. In particular, the 
exchange interactions become explicitly disorder depen- 
dent (Section IIA). Furthermore, by going beyond the 
Heisenberg limit, the operator for the Neel order param- 
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eter has to be corrected to take into account the charge 
mobility of the electrons in the Hubbard model 4 ^. The 
results presented below show that this correction is cru- 
cially important to obtain the correct 1/S behavior of the 
model. The consequent reduction in the amplitude of the 
staggered magnetization in the presence of local disorder 
is discussed in Section IIB. We then discuss in Section 
IIC the stability of the classical Neel ground state for 
finite disorder, when ring exchange is present. Section 
IID describes the spin wave method that we use. Sec- 
tion III gives an overview of the algorithmic procedure 
used to diagonalize the quadratic form of the disordered 
finite-lattice spin Hamiltonian. The numerical results are 
presented in Section IV, followed in Section V by a dis- 
cussion of the results and a perspective for future work. 
An appendix discusses the question of statistical uncer- 
tainties in the data presented in Section IV. 



in a perturbation expansion, to study the magnetic exci- 
tations and the staggered magnetization in the Hubbard 
mode l 41 i 46 . The method, introduced by Harris et a/4I 
and developed further by MacDonald et aZi 39 i 48 ' 49 , relies 
on the separation of the kinetic part T of the Hubbard 
Hamiltonian into three terms that respectively increase 
by one (If), keep constant (To) or decrease by one (T_i) 
the number of doubly occupied sites. Specifically, one 
writes: 

T = "* E <**A,c%* =T 1+ T + T_i (2.3) 



Ti 



-« E 



(2.4) 



II. SPIN HAMILTONIAN AND REAL SPACE 
LINEAR SPIN WAVE CALCULATION 

A. Spin Hamiltonian 

We begin with the Hubbard Hamiltonian, H^: 



H H = T + V 



(2.1) 



The first term is the kinetic energy term that destroys 
an electron of spin a at site j and creates one on the 
nearest neighbor site i. The second term is the on-site 
Coulomb energy U for two electrons with opposite spins 
to be on the same site i and where n, a 



the occupation operator at site i. A site i substituted 
by a non-magnetic cation has q = 0, otherwise = 1. 
In the following we use the notation to represent a 
summation over the L 2 sites of the square lattice and ^2 
for a sum over the N = Y^iO- ~ e «) undiluted sites. The 
number of magnetic sites and hence of mobile electrons, 
N, at half filling, is thus configuration dependent. The 
average concentration of vacancies is 1 — [eijdisordcr = x. 
Similarly, represents a sum over neighboring sites 
and a sum over neighboring occupied magnetic sites. 
Below, a summation index with angular brackets; (...) 
in > denotes an ordered sum, taking into account 
only unique pathways. 

The derivation of a spin Hamiltonian from a one-band 
Hubbard model can be performed through many differ- 
ent methods, leading to apparently different effective spin 
Hamiltonians. It is only recently that it has been showr* 4 ^ 
that all these Hamiltonians are equivalent, as they are 
related to each other through a unitary transformation. 
We have recently applied the canonical transformation 
method, which uses the ratio t/U as a small parameter 



1 E ('<' A.-'W.,-''...- 



r_i 



+ n i,(T C i, C r C j, ( j n j,o- 



t Ey £ie 3^ li ^ C i,(T C j,a n j,<' 



(2.5) 
(2.6) 



where a stands for up if a is down and for down if a is up 
and where = 1 — n^-j. This separation comes from 
multiplying the kinetic term T on the right by n ljff + 
hi tS = 1 and multiplying on the left by n,j^ + hj^ s = 1. 

Applying a unitary transformation e lS to leads to 
a spin-only Hamiltonian through the relation: 



_iS tt a —iS 
C DhC 



[iS,H H ] , [iS, [iS, Hn}] 



1! 



2! 



(2.7) 

We do not reproduce the derivation here, rather we refer 
the reader to Refs. for details of the form of 

S and Hh order by order in the development. Up to 
third order in the t/U expansion, we finally find for the 
effective spin Hamiltonian: 

«i,fe» 

+ J2 J3(Si-S m ) (2.8) 

«W)) 

+ ]T JcKs, ■ Sj)(s k ■ so 

<«,j>W> 

+(Si ■ St) (S k ■ Sj) - (Si ■ S k ) (Sj ■ SO} , 

where the site labels refer to the configuration shown 
in Fig. [2] The different coupling constants (Ji, J2, J3, J c ) 
arise as a result of the integration over all electronic paths 
allowed in the site-diluted Hubbard model. As a result, 
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i j m 



1 k 

FIG. 2: Labels for the different sites involved in the effective 
spin interactions and arising from at — U Hubbard model up 
to order t(t/U) 3 . ((i,k)}, (({i, m})) are nearest, second 

nearest and third nearest neighbors, respectively. (i,j,k,l) 
denotes the sites that belong to an elementary square plaque- 
tte. 



neighbor sites, these interactions are progressively elimi- 
nated as intermediate sites are diluted. That is, if both 
sites j and I are missing then J2(i,k) = 0. Hence, one can 
see that site dilution strongly affects the coupling con- 
stants as further neighbor exchange depends on the ex- 
istence of a nearest neighbor pathway between the sites. 
This would not be the case if the original Hubbard model 
included direct second or third nearest neighbor hopping 
parameters, t' and t", respectively^. We will return to 
this issue in the Conclusion section. However, in this pa- 
per we limit ourselves to nearest neighbor hopping only. 



they depend on the local site occupancy along the ex- 
change path. We find: 



Jo 



H 2 
u 3 



u 3 



(4e l£j +N lj ) 



jj^ (ei€je k + titit k - Nik, 



(2.9) 



(2.10) 



B. Neel order parameter 

Our objective is to calculate the ground state Neel or- 
der parameter for the original Hubbard model as a func- 
tion of site dilution, using a spin-only description. To 
do this, the staggered (spin density wave) magnetization 
operator, 



Aft 



nt), 



(2.14) 



J*=4 



J c = 80 



t 4 

IP' 

U 3 



(2.11) 



(2.12) 



where M^v is a plaquette index for bond \w and is equal 
to the number of plaquettes to which both sites \x and 
v belong. When there is no dilution, M^ v = 2 for all 
nearest neighbor bonds and 1 for second neighbor 
bonds across the diagonal of a plaquette. When 

one of the four bonds defining a plaquette is missing, 
the Nfn, for the three nearest neighbor bonds along the 
remaining edges of the plaquette are reduced from two to 
one. M^v for the next nearest neighbor bond across the 
diagonal of the plaquette is reduced from one to zero. For 
example, consider Fig. [5] where only the site j has been 
eliminated by dilution. The expression for the coupling 
constants becomes: 



J x {i,l) 



r i 4 



J*{i,k) = 4 



(2.13) 



U 3 



J 3 (i,m) = J c (i,j,k,l) = 0, 



which should be compared with J x = U 2 /U - 24t 4 /U 3 , 
J 2 = J 3 = 4t 4 /U 3 and J c = m A /U 3 for the undiluted 
lattice. The most important point here is that since 
the antiferromagnetic and frustrating J 2 and J3 exist 
solely via electronic hopping processes connecting nearest 



defined for the Hubbard model, must be canonically 
transformed before it can be exploited in a spin-only 

description. Here, Mh, M s and M s refer to opera- 
tors while M s and M s refer to their expectation val- 
ues. That is, within the effective theory Mh becomes 
M s = e lS Mue~ lS and the expectation value in the 
ground state is defined 



h(0|M h |0)h 
h(0|0) h 



3|M B |0) B 
,(010). 



(2.15) 



Here |0) H and |0) s 



|0)h are the ground state wave 



vectors in the original Hubbard and spin-only models. 
We have recently shownii that this is more than just an 
academic point. Rather, it has important consequences 
for the ground state magnetization as one moves into 
the intermediate coupling regime and, as we will show 
below, plays a significant quantitative role in the present 
site-diluted Hubbard model. As we apply the canonical 
transformation on M tt 41 i 46 we find for AL: 



M s = M H + i ( ! 
where 



1 



1 

2U~ 2 



T-iTt - T_iTi 



N n = [t u m h ], 

1 f_ 1 = [T_x,Mjj], 



N 



lf = [Tb, Mh]. 



(2.16) 

(2.17) 
(2.18) 
(2.19) 
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After some algebra, we can write this expression in terms 
of S = 1/2 spin operators 4 ^ as: 



AL 



N 
2t 



NU 2 



J2^j{s! -s*}(-iy. (2.20) 



(id) 



Recalling the standard definition for the staggered 
magnetization operator in a spin model, 



N 



we arrive at the principal result of Ref. [4jJ that 



(2.21) 



3 (0|A/ S |0) B 
3 (0|0) s 

s (0|M s |0) f 



s (0|0) s 



and 



(2.22) 



The difference is due to the appearance of new quan- 
tum fluctuations arising from the charge derealization 
over closed virtual loops of electronic hops, which is the 
origin of the second term in Eq. (|2.20p . These spin inde- 
pendent fluctuations, which appear to order t 2 /U 2 in the 
magnetization operator, are generated when the canoni- 
cal transformation is applied on the Hamiltonian to or- 
der t 4 /U 3 and are therefore not present in the (t/U — > 0) 
Heisenberg limit. We recently investigated the effects of 
these terms in the undiluted cas o 41 i 46 . Here, the disorder 
is manifest through the dilution variables e$ and, in this 
paper, we are interested in how the spin renormalization 
factor modifies the ground state magnetization upon site 
dilution. However, before doing so, we first return to a 
discussion of the ground state of the spin-only Hamilto- 
man Hs ■ 

Henceforth, for the sake of compactness, we shall omit 
the subscript "s" in M s and M s , understanding that all 
results presented below were obtained from calculations 
performed on a spin-only description of the low-energy 
sector of the half-filled Hubbard model. 



C. Classical ground state 

1. Ji interactions only 

The real space spin wave method that we use to deal 
with dilution requires, as the starting point, the knowl- 
edge of the classical ground state spin configuration. 
With nearest neighbor interactions only, the classical 
ground state configuration is, in the absence of dilution, 



the Neel staggered spin configuration. This long range 
ordered state results from the local minimization of the 
exchange interactions. Since we work with a concentra- 
tion of defects, or dilution, x, smaller than the perco- 
lation threshold x p , there exists a percolating cluster of 
magnetic sites with an exchange path connecting every 
pair of spins on the cluster. As a result, the classical 
ground state configuration for the percolating cluster is 
a connected Neel configuration, where every spin keeps 
the orientation it would have had without dilution (see 
Fig.©. 




FIG. 3: Diluted Neel configuration. The circle labels a miss- 
ing (diluted) site. 



2. Full Hamiltonian 

In the case of the effective spin-only Hamiltonian, ex- 
pressed in Eq. (|2.8p . the situation becomes more compli- 
cated. If the J2, J3 or J c interactions get too large, the 
system undergoes a phase transition to a new classical 
state that is not collincar. 

a. Non diluted case 

As can be read from Eqs. (|2.9|2.10|2.12j) . when there 
is no dilution, the coupling constants read: 



h 
■h 

J, 



= 4- 



24- 



t l 



= j 3 = 4-= 



t 4 

8 v 



u 3 



(2.23) 



For t/U = 1/8, a value similar to that reported for 
La2Cu04 and that we henceforth take in the present 
wor k 35 i 46 , the ratios between the different coupling con- 
stants are: 



r j_2 

Ji 
■h 

c 



0.0172, 
0.0172, 
0.0862. 



(2.24) 



For a model with nearest and next nearest couplings only, 
the J1/J2 model, the Neel state is stable for J2/J1 ^ 
0.5 50,51,52 p or tnc j 1 jj c model, the quantity J c = 
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J C S 2 is usually introduced, and as long as J c ^ J/2, 
the Neel state is stable 5 ^. Our parameters are far away 
from these critical values, and hence the classical ground 
state, without dilution, is Neel ordered. 
b. Diluted case 

One might have expected that the combination of frus- 
tration, brought about by J2 and J3 and site dilution 
would trigger an instability in favor of a local Villain 
canting of the spins^i, leading ultimately to a two di- 
mensional Heisenberg spin glass before x p is reached^. 



We note, however, that La2Cu04 is only approxi- 
mately described by the one-band Hubbard model with 
nearest neighbor hopping only. For instance we have 
recently shown that one can achieve a quantitative im- 
provement to the fitting of the spin wave excitation spec- 
trum measured by Coldea et al£^ by including direct fur- 
ther neighbor hopping constants t' and t" — . Such direct 
hops could change the above results, leading to canted 
classical ground state a 21 i 22 ' 23 i 24 ' 25 i 26 before the percola- 
tion threshold is reached (x c < x p ). 




FIG. 4: Particular dilution configuration. In this example, 
sites 2 and 4 are removed by dilution. 



However, as alluded to in the discussion below Eq. (12.131) . 
such locally Villain canted states do not occur in the 
model considered here, where all effective magnetic inter- 
actions derive from electronic processes involving nearest 
neighbor hopping. Hence, as we saw in the previous sec- 
tion, for the configuration of diluted sites shown in Fig. d] 
the second neighbor interaction, J13, between sites 1 and 
3 is destroyed by the dilution of sites 2 and 4. As a conse- 
quence, as long as the critical ratios for the J2/J1, J3/J1 
or J C J Ji for destroying two sublattice collinear Neel order 
arc not reached, there are no spins coupled by dominantly 
random frustrating interactions, J2, J3 or J c , as can be 
verified from studying Eqs. (|2.9l2.1QI2.12p . 

We therefore conclude that the classical ground state of 
7f s (4) in Eq. (EU) for t/U = 1/8 on the percolation clus- 
ter is a Neel configuration for all concentrations below 
the percolation threshold. From this, one can immedi- 
ately see the importance of the site percolation threshold 
in this problem: within the model considered, that is 
the site-diluted one-band Hubbard model of Eq. 12.21 the 
only accessible classical ground state is Neel ordered all 
the way to the percolation threshold x p . Hence, any re- 
duction in the range of stability of the Neel ground state 
is due uniquely to quantum fluctuations and is not due to 
(classical) random frustration effects. This conclusion is 
explicitly verified post factum within the real space spin 
wave calculation: any instability towards a non-collincar 
ground state would be detected as a negative eigenvalue 
of the Hessian matrix leading to complex eigenfrequen- 
cies. No such instabilities were detected in more than 
the ten thousand realizations of disorder considered in 
this work. 



D. Elementary excitations of a diluted spin 1/2 
system 

1. Method 

The introduction of site dilution destroys translational 
invariance, which excludes the use of Fourier space for 
calculating the spin-wave excitations. Hence, we closely 
follow the method introduced by Walker and Waldstedt^ 3 - 
to study excitations in Heisenberg spin glasses. Other re- 
cent studies of site-diluted 5" = 1/2 Heisenberg antiferro- 
magnets have followed a similar approac h 44 ! 45 . We first 
summarize this method for the simplest case of nearest 
neighbor exchange, J\ only, with Hamiltonian 



(2.25) 



As we know the classical ground state of the system, we 
can define for each site i, a unit vector pointing in 
the direction of the classical spin Si in this state. Note 
that in Eqs. (|2. 1412. 2012. 2"Tj) a unique global quantization 
axis in the lab frame, z , was used to define M s and 
M s . Henceforth, we label the spin components in terms 
of the projection of Si along the axis of a local right 
handed frame. We do so to keep with the original nota- 
tion of Ref. [43[ , from which we borrowed the method we 
use here. Let {xi,yi,n®} be an orthogonal triad of unit 
vectors and let pf and p~ be vectors defined by 



P 



Pi 



+ _ Xj + lyj 



V2 ' 
x t - iyi 



(2.26) 



V2 



We also introduce spin deviation (boson creation and an- 
nihilation operators), ai and a\, defined by 



S 4 -n° = S-a\ ai 



Si-pf = V2S 
S iP 7 = V2Sa\ 



25 



(2.27) 
(2.28) 



~2S~ 
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where the spin components are defined with respect to 
the local basis set, {a5j,yj,n°}. With Eq. 12.271 and the 
definition of pf , we can rewrite the Hamiltonian Eq. 12.251 
to order O(S) as: 



(id) 



1„. 

2 



(id) 



(2.29) 



+ n° • pfa] + n° ■ p t a l 

-S J ij [(Pi4 + Pi a i) ■ (Pj +a j + Pj a l) 

(id) 

- n° • n°(a\ai + aja,) . 



By making reference to the classical ground state, we 
introduce A^ defined by 



(2.30) 



Physically, A^ corresponds to the local staggered mean- 
field at site i originating from all the spins Sj to which 
Si is coupled. This change of variables makes the second 
term of the Hamiltonian vanish. We keep only the leading 
quantum correction to the classical term ^S 2 X^ij) n i ' 

n°, H2, quadratic in the {aj,^}: 



Ho = S 



. i 

\ J « ( p t a \ + P7 a i) ■ (Pt a ] + P. 



(id) 



The quantum-mechanical equations of motion are: 
da] 



(2.31) 



dt 



H 2 ,a. 



(2.32) 



.dai 

~ l ~d7 = ^ 2,a ^ ' 



which can be written: 



da; 



5 = -MEM + E p - 



(2.33) 



da 



where /', , = \,<),, - J,,\p, ' ■ p, ) and Q,, = J, , ; p, 1 • 
Pj + )- We use a vector representation for the operators 



and a\ ; that is a and a^ are iV-dimensional vectors whose 
components are and a\, respectively. As N, the total 
number of (occupied) magnetic sites, is configuration de- 
pendent, so are all vectors and matrices in the following 
discussion. We write 



d ( a 
eft V a* 



-P -Q 
Q* P* 



(2.34) 



where we refer to P and Q as the "interaction matrices" 
of order N x N. We can also write 



H 2 = S (a 1 a) H 



,t 



(2.35) 



where the Hamiltonian matrix H is defined in the 2N x 
2N phase space to be: 



H = 



-P -Q 
Q* P* 



(2.36) 



In order to diagonalize H, we perform a Bogoliubov 
transformation that introduces new boson operators d 
and d\ as follows: 



a = g*d + fd ] 
a< = f*d + gd^ 



(2.37) 



/ and g arc N x N matrices that must satisfy the boson 
commutation rules: 

g*f T -fg^ = 

9*9 T -fP = 1 , 

where f T is the transpose matrix of /. We can also write 
these relationships in matrix representation: 



EI& = I 



where 



E 



9* f 

r 9 



-1 
1 



(2.38) 



(2.39) 



are of dimension 2N x 2N . 

The aim of the Bogoliubov transformation is to diago- 
nalize Eq. (|2.34p . Consequently we require 



dt \ d f 



-ft 
ft 



(2.40) 



where ft is a diagonal matrix of eigcnfrcqucncies. Using 
(|2~34|) and (|2~39)) one obtains: 



If 9* f\fd 
dt\r 9 ) U f 



9* I 
f* 9 

-P -Q 
Q* P* 



-ft 
ft 



d 
rft 



9* f 

r 9 



d 



(2.41) 
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Hence, the equation we ultimately have to solve is: 

ED = HE, (2.42) 

where we have defined the complete matrix of eigenval- 
ues: 



D = 



-n o 
o n 



(2.43) 



With this method, we can calculate the zero point quan- 
tum spin fluctuations to order 1/5, and hence the expec- 
tation value for the spin on occupied site i: 



(Sf) = S - (4a,) = S 



\f„ 



(2.44) 



With the expectation value (S 1 ?) now defined in terms 
of |/iy| 2 , one can calculate the staggered magnetization, 
defined in either Eq. (f2T20|) for the finite t/U Hubbard 
model or Eq. (|2.21|) for the t/U — ► Heisenberg model. 
Formally speaking, in a thermodynamically large system, 
spins that reside on finite-size clusters and which are con- 
nected to the percolating cluster do not participate to the 
symmetry breaking nor do they contribute to the average 
bulk staggered magnetization. Hence, to capture that 
physics in the present problem, and to proceed numeri- 
cally, numerically, we first identify for a given realization 
of disorder, a percolating cluster of sites connected via 
nearest neighbor hopping. For each spin on the percolat- 
ing cluster, (Sf ) is determined from Eq. (|2.44[) . summed 
over, and normalized by N, the total number of sites for 
that realization of disorder (percolating and not), to give 
M 8 in Eq. (|2~2"0)) (henceforth denoted M). One then re- 
peats the calculation for many dilution configurations for 
a given x, performing a disorder average and obtaining 
both the averaged staggered magnetization on the perco- 
lating cluster, [M(x)] porc , or the bulk staggered magne- 
tization, [M(x)], averaged over all magnetic sites in the 
sample. We stress that, while the staggered magnetiza- 
tion on the percolating cluster, [M(a;)] per c is the most 
relevant quantity for the numerical study, it is the aver- 
age staggered magnetization over all Cu magnetic sites 
in the system, percolating and not, [M(a;)], which is ac- 
cessible to experiment, and which is displayed in Fig. [T] 
In the presence of interactions beyond J\ (i, j), the only 
change in the details of the above method occur in the 
matrix elements of P and Q. The form of these matrices, 
taking into account the second (J 2 ), third (J3) and ring 
(J c ) exchange interactions is discussed next. 



2. Calculation of the P and Q matrices 
a. Ji: first NN 

In this case the quadratic Hamiltonian reads: 



a]a l + OjCtj - a t aj - ajaj 



(id) 



(2.45) 



The P and Q interaction matrices then have the following 
form: 



/A, 



P 



\ 

Q = 



\ 



\ 



-Ji 



-Ji 



X L 2 J V J 

(2.46) 

where Xi is defined in Eq. (|2.30[> . and Q is a symmetric 
matrix = —Ji(i,j)- 

b. J 2 : Second NN 
We have: 

H 2 {J2) = S ^ J 2 (i,j) fa| a < + a ) a 3 ~ a \ a j - a ) a i) 1 

(2.47) 

which leads to the following additions to the P and Q 
matrices: 



( K 



(2) 



\ 



)(2) 



0/ 



(2. 



Ap is defined in a similar way to A,: 



(2.49) 



where ((j)) indicates a sum over the second neighbors of 
site i. 

c. J 3 : Third NN 



H 2 (J 3 ) = S ^ J 3(i,j) \ a\ai + ajaj - a\a.j - aja, j 
(((id))) 

(2.50) 

hence the expression for the P and Q matrices are mod- 
ified by: 



(3) 



p(3) 



Q 



(3) 



/o 



0/ 



(2.51) 
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with A| 3 ' defined by: 



«0»> 



°-^ (3) n°, 



(2.52) 



where (((?))) indicates the third neighbors of site i. 

d. J c : Ring exchange interaction 

To first order in l/S, the four spin terms appearing in 



J 



the Hamiltonian arc decoupled into bilinear products of 
a\aj. That is, to order l/S, the net effect of the ring 
exchange is to simply renormalize the J\ and J2 intcrac- 
tion o 35 ' 36 i 41 . The contribution of the ring exchange terms 
to the quadratic Hamiltonian is thus: 



H 2 {Jc) = -S 3 ^2 Jc(i,j,k,l) (a\ai + a^aj + a\a k + a\a^ + (a\a k + a\a t + a^ai + a\aj 

(i,j,k,l) 

- (a]aj + a]flj + a|a ; + a\a t + aja& + a\aj + a\ai + a\a k 

I 



where J c (i,j,k,l) = e^e/ce; J c . The elements of the Bo- 
goliubov transformation matrices g and / are thus mod- 
ified by the configuration dependent renormalization of 
the first and second neighbor exchanges. In zero dilution, 
Ji and J2 are renormalized to 35 ' 36 i 41 : 



— J\ — 25* J c — J 1 



J| ff — J2 — S 2 J C — J-2 



2 

Jr 



(2.53) 



III. ALGORITHMIC CONSIDERATIONS 

In order to obtain the quantum magnetization correc- 
tions in the disordered lattice, we have to solve the eigen- 
value problem described in Eq. (|2.43p . Results in the 
thermodynamic limit are estimated by doing a finite size 
scaling analysis for different system sizes. For each value 
of size and dilution, we generate many realizations of dis- 
order after which we perform successively the disorder 
average and the finite size scaling to the thermodynamic 
limit. Our algorithm is organized as follows for each value 
of the system size and dilution: 

• Generation of the diluted lattice and computational 
identification of the percolating cluster. 

• Calculation of the P and Q matrices (Eq. (|2.34D ~). 

• Diagonalization of the matrix using Lapack rou- 
tines. 

For a system of linear size L and dilution concentra- 
tion x, for each site, we generate a random number r 
between and 1. The site is considered as removed if 
r (1 — x). For each realization of disorder, for which 
the number of sites N(L 2 ,x) is different, we first con- 
struct the percolating cluster. To do this, the undiluted 



sites are labeled from 1 to N. Starting from site 1, with 
coordinates (a, /?), we verify if the neighbors (a ± l,/3), 
(a,/? ± 1) are occupied. If yes the label of the site is 
changed to 1 . Moving to one of these sites the procedure 
is repeated. If the cluster 1 terminates, the next cluster 
takes the number of the first occupied site encountered. 
Once all sites have been visited, the procedure is repeated 
taking an arbitrary starting point. If neighboring sites 
are occupied the indices of the two sites take the lowest 
of the two values. The procedure is repeated until no 
further evolution occurs. For the biggest cluster we then 
check for the existence of percolating pathways along the 
x and y direction. If a percolating cluster exists, then 
the matrix H (|2.36p is constructed. 

The diagonalization of H is performed using a fortran 
77 Lapack double precision set of routines: 

• DGEHD2 computes Hessenberg reduction of the 
H matrix. 

• DORGHR and DHSEQR lead to the Shur fac- 
torization. 

• DTREVC gives the eigenvectors of the H matrix. 

From the results of the Lapack routines we first con- 
struct a matrix of eigenvectors E of H . We order the 
columns of E so that the first column is an eigenvector 
corresponding to the lowest eigenfrequency, and the last 
column is an eigenvector corresponding to the highest 
eigenfrequency of H. 

The matrix E is thus defined up to the subspaces of 
degenerate eigenvectors and the matrix D in Eq. (|2.43[) 
is the diagonal matrix of its cigenfrcquencies: 



ED E~ 



H 



(3.1) 



However, knowledge of D and E does not completely 
determine the problem. In order to establish the elements 
of the Bogoliubov transformation we must construct from 
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E the matrix E that satisfies both the relation (|2.38[) . 
coming from the boson commutation relations and the 
eigenvalue Eq. (j2T43|) . That is: 

El3 = I and EDE- 1 = H (3.2) 

We find E through the application of a transformation: 

E = Edb, (3.3) 

where db is a block-diagonal matrix. Using the commu- 
tation relation (|2.38| one finds 

dbldb^iAiy 1 (3.4) 

with 

M = E ] IE. (3.5) 

M is a Hcrmitian matrix obtained from the Lapack rou- 
tines. It is block diagonal, with blocks Mi of size Pi x p il 
corresponding to a subspace of degenerate eigenvalues of 
H, of dimension pi. The transformation matrix db is 
therefore also block diagonal, with corresponding blocks 
Ai. If the matrix, Ei, represents the subspace of eigen- 
vectors, of dimension 2N x pi, the transformation gives 
Ei = EiAi. 

To find Aj, we need to solve 

±AU] = Mr\ (3.6) 

where the sign ± depends on which sector of the eigen- 
value matrix D in Eq. (|2.43[) the subspace belongs. The 
blocks Mi are first inverted and then diagonalized 

Mf 1 = Ki DiKr 1 . (3.7) 

The matrix Di contains either positive or negative eigen- 
values, depending on the sign of the eigenvalue of the 
subspace of H. From this we find the diagonal matrix 
\/±Di, where the sign ± is chosen so that the square 
root is defined. 

Ai is finally found from: 

Ai = Ki^±DiKr\ (3.8) 

In this block diagonal procedure the subspace corre- 
sponding to the Goldstone modes is explicitly excluded. 
All other operations are then mathematically well de- 
fined^ and the Bogoliubov transformation is completely 
determined. A further summary of the calculation pro- 
cedure can be found in Appendix [A] 

IV. RESULTS 

We have calculated the quantum fluctuations of the 
magnetic ground states for a diluted spin system for two 
situations: 



1. t/U — > 0, Hcisenberg limit. The result are com- 
pared with those obtained by Mucciolo et alM- who 
used a similar linear spin wave calculation. The 
motivation here is to validate the two sets of re- 
sults against each other and to quantify finite size 
effects and statistical errors. In this limit the addi- 
tional terms in the magnetization operator (|2.16p 
leading to the inequality (|2.2ip are zero. This cal- 
ibration allows us to confirm that there is indeed 
a discrepancy between the results from 1/5 cal- 
culations^ and those from quantum Monte Carlo 
simulations^ in the Heisenberg limit. 

2. t/U = 1/8. This t/U value is close to the one 
found for La 2 Cu04 by Coldea et al^- (t/U ~ 
1/7.35)^. By using this value, we can begin an 
investigation of the effects of further neighbor and 
ring exchange interactions in the experimentally 
relevant situation of Cu substitution by Mg and 
Zn in La2Cui_ x (Mg/Zn)j;04. 

A. t/U = 0: diluted Heisenberg model 

The presence of dilution introduces statistical fluctu- 
ations in the ground state magnetization due both to 
configurational variations for fixed number of magnetic 
sites and to the variation in the number of magnetic sites 
from one configuration to another. To combat this, we 
perform an average over a number of disorder configura- 
tions, Ao, that increases with the dilution. We chose Ao 
to be the integer closest to 5000 times x, where x is the 
concentration of missing magnetic sites. Averaging over 
the disorder we define the average number of sites for a 
given concentration: 

R = jf T, N ( L2 > x ) = (!-^ (4- 1 ) 

where N(L 2 , x) is the number of sites for a system of size 
L and concentration x for a specific disorder realization. 
System sizes were studied from L 2 = 10 x 10 to L 2 = 
34 x 34. A detailed discussion of the various contributions 
to the statistical errors can be found in Appendix [Bl 

The ground state magnetization is estimated by ex- 
trapolating the finite size results to the thermodynamic 
limit. In order to do this we proceed in two steps: 

• Firstly we determine the staggered magnetization 
for the sites on the percolating cluster for each 
realization of disorder; from which the staggered 
magnetization averaged over all magnetic sites for 
a specific realization of disorder is obtained. We 
then make a disorder average over many real- 
izations, calculating both the disorder averaged 
staggered magnetization on the percolating clus- 
ter, [M (x, L)]p C rc, and the experimentally relevant 
disorder averaged bulk staggered magnetization, 
[M(x, L)]. The errors on these measures are es- 
timated as explained in Appendix [Bj 
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• This process is repeated for different system sizes 
for a given x and the results are extrapolated to the 
thermodynamic limit by making a least squares fit 



[M(x,L)} pelc = [M(x)] pcrc + - + — . (4.2) 
The same procedure is used for [M(x, L)]. 



As an example, we show results for x = 0.2 in Fig. 
[S] where we plot the magnetization [M (x, L)] porc against 
1/L. As one can see, the statistical noise on the data is 
small and is consistent with the size of the error bars es- 
timated in Appendix [Bl The magnetization extrapolates 
linearly to the thermodynamic limit in 1/L to an excel- 
lent approximation, allowing a high precision estimate 
for [M] pcrc : 



[M(x = 0.2)]j 



0.236 ±0.001 



(4.3) 



We note that between L = 10 and the biggest system size 
studied, L = 34, [M] pe rc varies by over 30%. This sub- 
stantial variation confirms the need for such a finite size 
scaling procedure here. Results for different values of x 
are shown in Fig. [6] For the system sizes studied, the size 
dependence is very nearly linear in 1 / L for all x. One can 
also notice that the slope, a, is almost independent of x 
until the percolation threshold, x p = 0.41 is approached, 
at which point it increases with finite size effects becom- 
ing progressively more important. This evolution is not 
inconsistent with the critical nature of the percolation 
threshold and the question as to whether [M] porc , de- 
termined via the 1/S method, goes continuously to zero 
or jumps discontinuously to zero at x p is an intriguing 
one. On the other hand, it is found from quantum Monte 
Carlo simulations that [A/] pcrc has a discontinuous jump 
at x = x p — . However, this question is not the main 
focus of the paper and to do it justice would require a 
more extensive and dedicated study near x p . Here we 
simply remark that [M] perc extrapolates to small values 
for concentrations less than, but near x p . 

Collecting these results, we show the staggered mag- 
netization for the ground state of the site-diluted Hciscn- 
berg model, as a function of dilution in Fig. [JJ [M] pcrc 
goes smoothly from the known value for the undiluted 
case in the 1/S approximations^^, [A/] pcrc « 0.31, 
to zero for x very close to the site dilution percolation 
threshold, x p . 

In Fig. [51 we compare our result with those obtained 
by Mucciolo et a/4i for the same model. The data are 
normalized by the value [M] perc (x = 0) = M(0). There 
is extremely good quantitative agreement between our re- 
sults and theirs, providing strong evidence that the two 
methods give correct results for the 1/S* method consid- 
ered. 

It is important here to make a comparison between our 
results and those from quantum Monte Carlo (QMC), 
which is in principle exact, apart from numerical error. 
Such a comparison is made in Fig. [5] where we show 
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FIG. 5: Staggered magnetization for the t/U — > Heisenberg 
model for x = 20% and L G [10, 34]. 
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FIG. 6: Evolution of [M] pcrc for the t/U — > Heisenberg model 
with L for different values of dilution — from top to bottom: 
x = 0%, x = 6%, x = 16%, x = 26% and x = 36%. 



unnormalized data for the magnetic moment on the per- 
colation cluster from our calculation, compared with the 
QMC data of Ref.[29|. For zero dilution, the methods give 
very similar results. This is expected as it is known that 
1/S 2 contributions to the quantum fluctuations in this 
case are identically zero 58 i 60 ' 61 , meaning that the differ- 
ence between spin wave and QMC comes, to leading or- 
der, from 1/S 3 contributions, which one might expect to 
be small. Moving away from zero dilution, the difference 
between the two sets of results increases in a monotonic 
way, with the moment from the QMC consistently larger 
than that determined from the 1/5 spin wave calcula- 
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FIG. 7: Staggered magnetization for the t/U — > Heisenberg 
model on the percolating cluster extrapolated to the thermo- 
dynamic limit. The solid line is a guide to the eye. 
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FIG. 9: Staggered magnetization on the percolating cluster 
for the site-diluted Heisenberg model: comparison between 
QMC (from Ref. [2^]) and 1/5 spin wave results. The solid 
and dashed lines are guide to the eye. 
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FIG. 8: Comparison between our data for the ground state 
(bulk averaged) staggered magnetization for the Heisenberg 
model, normalized by the value at zero dilution, and the data 
of Mucciollo et aZ.— The dashed-dotted line is a guide to the 
eye. 



tion. Hence the comparison explicitly illustrates that the 
1/5" method over estimates the importance of the quan- 
tum fluctuations in the presence of disorder. In order 
to understand this quantitative difference, one should in- 
vestigate the effects of magnon-magnon interactions and 
Berry phase terms, which we do not attempt here. 

The above limitations should be taken into consid- 
eration when comparing data from spin wave calcula- 
tions with the experiment. In Fig. [TU] we add our re- 



sults to those previously shown in Fig. 1 for compar- 
ison with the experimental neutron scattering data on 
La2Cui_ 2; (Mg/Zn) 2 ;0^ and quantum Monte Carlo sim- 
ulations^. The figure allows us to confirm the conclu- 
sion, already made in Ref. 0, that the extremely good 
agreement between experiment and the spin wave calcu- 
lation is rather fortuitous: if the Heisenberg Hamiltonian 
was an adequate starting point to describe the experi- 
mental data, the "exact" QMC results would be in bet- 
ter agreement with the experimental data than the l/S 
spin wave data are. As can be seen from the figure, the 
reverse is true; while the QMC data is consistently above 
the experimental curve, the spin wave data lies very close 
to it. Hence, although the Heisenberg Hamiltonian is 
clearly a good starting point for acquiring an acceptable 
qualitative description of Mg and Zn doped La2Cu04 it 
appears, on the basis of the results shown in Fig. 1101 
to be inadequate for a really quantitative description. 
Further, we remark that the experimental data are pre- 
sented such that they are normalized by [M(x = 0)] pcrc . 
While the undiluted moment is [M(x = 0)] pcrc w 0.31 
from QMC and spin wave calculations, recent estimates 
by Lee et alr& place the experimental moment at about 
0.25. Hence, removing the absolute scale of the magnetic 
moment improves the impression of a good agreement of 
the experiment with the QMC results for the site-diluted 
Heisenberg model for small x. When plotted on an abso- 
lute scale the agreement between experiment and theory 
would be less convincing. This is an important point 
for the present paper as we continue to work within the 
linear spin wave approximation and cannot expect to ac- 
count for the contributions beyond l/S linear spin waves 
which, following the QMC approach, appear to be im- 
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portant for the dilution problem. That said, by making 
improvements to the starting spin-only description of the 
Hubbard Hamiltonian, through higher order terms in the 
canonical transformation, we can expect to improve the 
comparison with experiment on an absolute [M(x)] scale. 

With this in mind we have extended our calculations 
to order t 4 /U 3 , which allows us to include second and 
third neighbor exchange as well as ring exchange around 
an elementary plaquette, and also to include quantum 
fluctuations from charge dclocalization in the underlying 
Hubbard model. 




Dilution x (%) 

FIG. 10: Ground state magnetization as a function of dilution 
for Mg and Zn doped La 2 Cu0 4 : La 2 Cui_ ;c (Mg/Zn) ;c 04^, 
for quantum Monte Carlo— for the site-diluted square lat- 
tice Heisenberg antiferromagnet (SLHAF). The dashed line 
is a guide to the eye parameterized by \M(x)]/M (0) = (1 — 
x/x P Y eii . The figure is reproduced from Ref. [30]. Added to 
this figure is our data and that of Ref. [3 obtained from the 
numerical 1/5 spin wave analysis of the site-diluted SLHAF. 



B. t/U = 1/8: on the role of the ring exchange 
interaction 

When interactions beyond nearest neighbor exchange 
are taken into account, two effects have to be consid- 
ered. First the transverse spin fluctuations are modified 
by the inclusion of the new interactions since these af- 
fect the magnon excitation spectru m 35 ' 46 . Secondly, the 
charge dclocalization induces a further quantum fluctu- 
ation term over and above those from transverse spin 

fluctuations. This is the difference between M s and M s 
in Eq. (|2 . 20[) and which leads to renormalization of the 
staggered magnetization in a way that depends on dilu- 
tion. In this section, we treat these two effects separately 
to quantify their respective importance for t/U = 1/8. 



1. In the absence of charge mobility renormalization 

Finite size results: 

The first point we wish to illustrate here is the impor- 
tance of the modification of the exchange pathways in the 
diluted system. We have argued in Section pi C| that 
dilution does not introduce random frustration at the 
classical level, even in the presence of further neighbor 
spin interactions, if these interactions are derived from 
the Hubbard model with nearest neighbor hops only. In 
this case, such a longer range interaction depends on the 
presence of a nearest neighbor exchange pathway. Hence, 
we do not expect long range interactions to have a desta- 
bilizing effect on classical Neel order on the percolating 
cluster. This can be seen indirectly by comparing the 
finite size scaling of our effective spin-only model with 
that of a more phenomcnological model. In the latter, 
which we refer to as the "p-model" , the further neighbor 
interactions have full strength, independently of the ex- 
istence of a nearest neighbor exchange path created by 
the electronic hopping processes, so that they exist even 
if the pathway is severed by a non-magnetic defect (i.e. 
diluted site). In the p-model, the bilinear exchange in- 
teractions Ji and J3 are taken to be J2 = 4i 4 /U 3 €iCj and 
J3 = 4t 2 /U 3 Ciek while J c is kept to have the same site oc- 
cupancy dependence as in Eq. (|2.12p . In Fig. [Til we show 
results for the size dependence of the staggered magneti- 
zation [iWjperc a s a function of concentration x of diluted 
sites, for the effective spin-only Hs in Eq. (|2.9p with 
{Ji, J2, J3, Jc} coupling constants as given in Eqs. (|2.9p 
to (|2.12p . The magnetization is a monotonic function of L 
for all values of x. This should be compared with Fig. [12] 
where we show similar data for the p-model. For large 
dilution, the statistics are much worse and the magneti- 
zation considerably lower than in the first case. This indi- 
cates the build up of random frustrated plaquettes that 
eventually destroy the Neel order before x p is reached, 
even at the classical leve l 21 ' 22 ^ 3 -^^ 5 -^ 6 -. 
Thermodynamic limit 

For the effective spin-only Hamiltonian, results are ex- 
trapolated to the thermodynamic limit, using the proce- 
dure described in the previous section and Eq. (|4.2p . In 
Fig. [TjJ]we show the ground state magnetization, [M] perc , 
compared with the previously shown results from Fig. [7] 
for the (t/U — >, Ji only) Heisenberg model. 

The first thing to notice is that there is very little dif- 
ference with the Heisenberg model! The second is that 
the small difference that is present is towards a higher 
ground state magnetization, with the maximum change 
occurring at x = 0%. This is because the ring exchange 
terms decrease the transverse spin fluctuations 41 "^. As 
discussed in Refs. [4lU4r| . this increase in the magnetic 

moment occurs because the ring exchange terms in 
decouple in the 1/5 expansion into effective ferromag- 
netic second neighbor two-body exchange terms which 
further stabilizes the two sublattice Neel order by reduc- 
ing the transverse spin fluctuation s 41 ' 46 (see Eq. (|2.53p h 
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FIG. 11: Size dependence of the staggered magnetization 
of the effective spin-only Hamiltonian ■ Further neigh- 
bor interactions are generated through the existence of near- 
est neighbor electronic hopping pathways. Charge mobility 
renormalization effects are not included here. 



FIG. 13: Dilution dependence of the staggered magnetiza- 
tion of the effective spin-only Hamiltonian obtained from the 
Hubbard model - No charge renormalization effect. The inset 
shows a blow-up for x < 20%. 
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FIG. 12: Size dependence of the staggered magnetization of 
a "p-model" Hamiltonian where further neighbor interactions 
are not explicitly dependent on the neighbor exchange hop- 
ping pathways. 



For x greater than about 12% dilution, this stabilization 
effect is largely destroyed and the two curves merge up 
to the percolation threshold. This is explained by the 
fact that, as these interactions involve more than two 
sites, they are more sensitive to dilution than the nearest 
neighbor terms, and their effect becomes negligible long 
before the percolation threshold is reached. The effects at 
high dilution would be very different for the p-model (see 



Fig. 1 1 2[) with frustrating further neighbor interactions 
that are independent of the presence of nearest neighbor 
exchange pathways. However, in the context of compari- 
son with experimental results on La2Cui_ x (Mg/Zn) a: 04, 
such terms should only appear through electronic hop- 
ping over further neighbors in the Hubbard model. We 
have recently considered this problem in the absence of 
dilution 4 *, but extending this work to include dilution is 
beyond the scope of the present study. 



2. Finite charge mobility renormalization 

The charge mobility, or electron delocalization ef- 
fect, leads to a decrease in the magnetizatio n 41 i 46 (see 
Eq. (|2.20[) V However, the delocalization is also condi- 
tioned by the allowed nearest neighbor electronic hop- 
ping pathways and is consequently also dilution depen- 
dent. The finite size scaling of the magnetization, as de- 
scribed by Eq. (|4.2[) . is not changed qualitatively by this 
renormalization (not shown here) and the results are ex- 
trapolated to the thermodynamic limit, using Eq. (|4.2[) . 
As shown in Fig. Q3J there is a significant decrease in 
the magnetization compared with [M] porc without charge 
mobility renormalization, or with the Hciscnberg model 
(see Fig. [13]) . This difference is again reduced as the dilu- 
tion increases. For x = 0% the decrease is of the order of 
14%, whereas it goes down to 9.5% for x = 30% and goes 
towards zero at the percolation threshold. We conclude 
therefore that the charge delocalization term is a major 
contribution to the corrections found by extending, to or- 
der t 4 /U 3 , the canonical transformation of the Hubbard 
model into an effective spin Hamiltonian. It is explicitly 
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a property of the Hubbard model and is not present in a 
phenomenological spin-only model. It is therefore clear 
that care must be taken when using such phenomeno- 
logical spin-only models without directly considering the 
mobility of the underlying system of electrons when aim- 
ing at obtaining a quantitative description and compari- 
son between experiment and a microscopic theory. This 
is the main result of this paper. 
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FIG. 14: Evolution of the staggered magnetization with dilu- 
tion for the effective spin-only Hamiltonian with charge mo- 
bility effects included. The open circles are the data for the 
Heisenberg (t/U —> 0) model. The diamonds and squares 
are respectively the data for the spin-only representation of 
the t/U = 1/8 site-diluted Hubbard model with and without 
finite charge mobility renormalization included. 



proposed here can describe many of the magnetic features 
of La 2 Cui_ :c (Mg/Zn) ;c 04. 

In Fig. [TOJ we re-plot the data of Fig. [H] normalized 
by the ground state order parameter at zero dilution, 
M (x = 0) . The two data sets for the Heisenberg model 
and that for the spin model with further neighbor 
interactions with the effects of electron delocalization via 
virtual hopping included, lie on top of each other. Hence, 
the improvements brought in by developing the effective 
spin description of the Hubbard model up to four virtual 
hopping terms, including ring exchange are not evident 
when the data are normalized in this way. The data sets 
would therefore continue to give the same favorable, but 
fortuitous comparison with the experimental data of Vajk 
et al^L, as seen in Fig. [TO] and discussed in Ref. [44| . 
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FIG. 15: Data from Fig. 1141 normalized by the ground state 
order parameter at zero dilution, M(x — 0). 



C. Experimental considerations 

Figure [TJ] illustrates our main result concerning the 
comparison with experiment: inclusion of the charge mo- 
bility renormalization factor shifts the scale of magnetiza- 
tion downwards over the whole dilution range. Compar- 
ing results for zero dilution; for the Heisenberg model, 
the ground state moment is [M s ] 0.31, while experi- 
ment yields [M s ] oxp w 0.25 — . Hence a comparison of 
[M(a;)] data not normalized by M (x = 0) will show the 
Heisenberg model, either from spin wave, or from QMC 
to be above those from the dilution experiments. Includ- 
ing hopping processes to order t 4 /U 3 for t/U = 1/8, a fair 
estimate for La 2 Cu0 4 , one finds^^ [M B ] « 0.27. This 
is still above the experimental value, but it is clear that, 
taken altogether, the extra corrections arising from both 
transverse spin fluctuations and finite electron mobility 
away from the t/U — ► Heisenberg limit, have scaled 
the magnetization in the right direction. This is an im- 
portant result indicating that the perturbative methods 



V. CONCLUSIONS AND PERSPECTIVES 
A. Conclusions 

In this paper we have investigated the problem of site 
dilution in systems described by the half-filled one-band 
Hubbard model. We have extended the canonical trans- 
formation technique to calculate an effective spin Hamil- 
tonian up to order t 4 /U 3 , for magnetic site dilution x. 
We use a real space spin wave technique, linear in (1/S) 
to calculate the dilution dependence of quantum fluctu- 
ations on the staggered magnetization. Specifically, we 
considered two problems. We first studied the Heisen- 
berg t/U — ► limit, comparing our results with those 
from quantum Monte Carlo (QMC) studies^ on the same 
model. We confirm, to a high degree of accuracy, previ- 
ous results from Ref. 0, using a similar 1/S technique. 
Hence, our results also confirm a systematic deviation 
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between the QMC — and the spin wave results for fi- 
nite dilution. This difference, which is small for zero 
dilution, illustrates the dilution dependent generation 
of magnon-magnon interactions and Berry phase terms, 
both of which are neglected in the l/S spin wave calcula- 
tion. By comparing the QMC results and the l/S results, 
one concludes that these two effects work to stabilize the 
semi-classical two sublattice Neel order rather than to 
drive the system into an exotic quantum phase. Hence, 
while the spin wave technique predicts the magnetic mo- 
ment on the percolating cluster going to zero at, or very 
close to the percolation threshold, QMC simulations find 
a renormalized classical result, with the moment on the 
percolating cluster taking a discontinuous jump to zero 
at the percolation threshold. 

The second and main objective of our work was to 
investigate how corrections to a site-diluted spin-only 
Hamiltonian, originating from a site-diluted one-band 
Hubbard model affects the dilution dependence of the 
Neel order parameter (staggered magnetization [M(x)]). 
Underlying this question, was the goal of obtaining some 
information on the capacity of the one-band Hubbard 
model to describe data from experiments on site-diluted 
La 2 Cu0 4 , La 2 Cui_ 3: (Mg/Zn) a: 04^. 

Within the one-band Hubbard model, the best esti- 
mates from fitting to magnon excitation spectr a 35 i 46 give 
t/U ~ 1/8, placing the system away from the Heiscnbcrg 
limit and into the intermediate coupling regime, where 
higher order electron correlations need to be taken into 
account. Integrating out the kinetic degrees of freedom 
via a canonical transformation, implemented to order 
t 4 /U 3 , introduces second and third neighbor interactions 
into the resulting spin-only Hamiltonian, as well as ring 
exchange term from electronic pathways around a closed 
square plaquette. Calculation of the magnetization op- 
erator for the original Hubbard model, introduces at this 
order a charge mobility term that renormalizes the mag- 
netization below that obtained by considering the (triv- 
ial) definition of the staggered magnetization of a spin- 
only Hamiltonian given by Eq. (|2.21[) . Including all these 
effects, and staying within the linear (l/S) spin wave 
approximation, we find a reduced estimate for the mo- 
ment at zero dilution, w 0.27/ib compared with 0.31/ib 
for the Heisenberg model — . As further neighbor and 
ring exchange interactions mediated by four electronic 
hops require unbroken exchange pathways over length 
scales greater than nearest neighbor, their effects disap- 
pear well before the percolation threshold is reached. The 
net result is therefore that the evolution of the ground 
state moment as a function of dilution is qualitatively 
similar to that for the Heisenberg model, disappearing 
at the percolation threshold in the same way, but with 
the absolute scale renormalized downwards by 10% to 
15%. While the 1/5 method is subject to the limita- 
tions described above, our results clearly illustrate that 
for an ultimate detailed and quantitative understanding 
of the role of site dilution in a correlated electron sys- 
tem, such as La2Cui_ 2; (Mg/Zn) a ;04, the charge mobility 



effects must be taken into account. Such a description is 
beyond a spin-only model, decoupled from an electronic 
model describing the behavior of the strongly correlated 
electrons. 



B. Perspectives 

To expand on the work presented in this paper, it 
would be interesting to carry out further theoretical and 
numerical studies using a common calculation scheme for 
both the site-diluted Hubbard model expressed in the 
framework of a spin-only Hamiltonian with ring exchange 
and the site-diluted Heisenberg (t/U — > 0) model. How- 
ever, in the absence of a solution to the sign problem 
for frustrated quantum spin systems, quantitative results 
for the generalized dilution problem from large quantum 
Monte Carlo simulations remain inaccessible. 

Angle resolved photo emission spectroscopic (ARPES) 
experiments as well as ab-initio calculations on a number 
of copper oxide materials provide strong evidence that an 
effective one-band Hubbard model description of these 
systems must include direct hopping parameters t and t" 
to second and third nearest neighbor sites. Furthermore, 
such experiments and calculations indicate that these pa- 
rameters are not significantly smaller than the nearest 
neighbor hopping t, with t'/t ~ —0.3 and t"/t ~ 0.15. 
We have recently included direct hopping parameters t' 
and t" in a derivation of a spin-only Hamiltonian repre- 
sentation of the half-filled t—t' — t"—U Hubbard model^. 
As a result of these sizeable energy scales, our analysis 
of magnon excitation spectra in La2Cu04 reveal that the 
contributions from these parameters are of similar mag- 
nitude to the four hop (order l/U 4 ) processes for nearest 
neighbor hopping, which give rise to the ring exchange 
interactions studied in Ref. [H| and in the present paper 
(last term in ff s (4) of Eq. (3]5J>). 

A key result of Ref. [46[ is that the ground state stag- 
gered moment, approximately 0.235, is further reduced 
from the value found for the t—U Hubbard model, ~ 0.27, 
using the t and U values of Coldea et al£^. This value 
is closer to, but undershoots the experimental estimate 
of 0.25 Although this progression lies within the ex- 
perimental uncertainty, the detailed analysis of Ref. [46| 
suggests that the t — t' — t" — U Hubbard model is a much 
improved starting point for a quantitative description of 
the magnetic properties of La2Cu04. This conclusion is 
in accordance with ARPES studies and ab-initio calcu- 
lations on various cuprates. 

A natural extension of the work presented in this pa- 
per would be to investigate the role of t' and t" in the 
site dilution problem. In this model a large number of 
new ring exchange terms are generated and the further 
neighbor hopping terms allow for connected pathways for 
dilution concentrations above the nearest neighbor per- 
colation threshold. It seems likely that these extra terms 
would change the shape of the [M s (a;)] vs x curve, spe- 
cially for x close to x p , and hence change the qualitative 
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aspect of the results even when the magnetization scale 
is factorized out of the problem, as in Fig. [T5] 

The realization that t' and t" are important energy 
scales in a Hubbard model description of La2CuC>4 leads 
to an interesting experimental puzzle when considering 
the substitution of Cu 2+ by non-magnetic Zn 2+ and/or 
Mg 2+ . As discussed earlier in this paper, and as illus- 
trated by the convergence of the the results of Fig. [151 for 
the Hciscnbcrg model (Ji only) and the t/U = 1/8 Hub- 
bard model (Ji, J2, J3 and J c ), the dilution-dependence 
of the electronic hopping pathways leads to a crossover 
concentration x* (x* ~ 15% for t/U =1/8) above which 
the influence of the J2 — J3 — J c terms of order l/U 3 
has essentially vanished. However, the presence of direct 
t! and t" hoppings leads to additional frustrating second 
and third nearest neighbor exchange with J' 2 « 4(t') 2 /U 
and J3 w 4(t") 2 /U, respectively. Unlike the J 2 and J 3 
interactions generated by fourth order hopping processes, 
J' 2 and J3 do not depend on the presence of nearest neigh- 
bor pathways and hence are unaffected by the dilution 
(see discussion in Section HVB H and the one accompany- 
ing Fig. [T2|) . As there is now frustration which is inde- 
pendent of the existence of nearest neighbor pathways, 
one would expect that upon dilution there would be a 
proliferation of Villain canted state o 21 i 22 ' 23 i 24 ' 25 i 26 as the 
concentration of impurities approaches the percolation 
threshold x p . This could ultimately lead to a Hciscnbcrg 
spin glass phase for a dilution concentration x < x p . In 
this context, it is perhaps surprising that experiments 
find sharp (resolution limited) magnetic Bragg peaks in 
La2Cui_ 2: (Mg/Zn) x 04 all the way to x = Xp^. It would 
certainly be interesting to revisit this question and study 
in more detail the possibility of a spin glass phase devel- 
oping in La 2 Cui_ 2: (Mg/Zn) 2 .04 close to the percolation 
threshold. We note further that in the region close to the 
percolation threshold there is the possibility of a freezing 
transition of the transverse spin components only. Such 
a transition could be observable in nuclear quadrupolar 
resonance (NQR) or muon spin relaxation (muSR) exper- 
iments^ as were done sometime ago on L^Cui-^Zn^C^ 
— . However, in those early experiments 2 ^, it now seems 
likely that the then detected transverse spin freezing was 
driven by doped holes introduced by an imperfect control 
of the oxygen stoichiometry in I^Cui-^Mg/Zn^C^ 
30 i 31 . It would be interesting to repeat such NQR and 
muSR experiments on La 2 Cui_ a; (Mg/Zn) 2 ,04 samples of 
the same quality as those used in neutron scattering ex- 
periments of Ref. [3(j ■ 

Another effect that could be relevant for 
La2Cui_ a: (Mg/Zn) 2 ;04 is the local distortion of the 
lattice due to the small difference in the ionic radius 
between Cu 2+ and Mg 2+ or Zn 2+ This difference 

could lead to a local modification of the hopping 
parameter t in the neighborhood of a site where a 
Cu 2+ ion is replaced by a nonmagnetic ion (see Fig. 2 
in Ref. [H[). Such disorder- induced variations of the 
hopping parameters could then contribute to explain 
the difference between the experimental data and QMC 



data in Fig. [T] The importance of local distortions could 
perhaps be provided by local probe experiments such as 
muSR, NMR or NQR. This problem may also be con- 
sidered as a precursor to the study of disorder-induced 
static magnetism in cuprate superconductors^. In this 
case the inclusion of mobile holes makes it much more 
complicated, but the study of the diamagnctic dilution 
problem in La 2 Cui_ 2; (Mg/Zn) x 04 maintained at half 
filling could provide a useful framework on which to 
build. 

In conclusion, we have explored in this work the prob- 
lem of the evolution of the magnetic order in a spin- 
only representation of a site-diluted one-band Hubbard 
expressed in terms of a spin-only Hamiltonian, taking 
into account up to four hop processes. For a finite ra- 
tio of hopping constant to on-site Coulomb energy, t/U, 
the resulting spin Hamiltonian differs from the simpler 
site-diluted S =1/2 Heisenberg model, containing effec- 
tive exchange coupling beyond nearest neighbor as well 
as ring exchange interactions. The long range exchange 
interactions, the ring exchange and the renormalization 
of the nearest neighbor exchange depend specifically on 
the local random hopping pathways that remain unin- 
terrupted by the missing (diluted) sites. We hope that 
this study can motivate further analytical and numeri- 
cal studies of the site-diluted one-band Hubbard model 
as well as new experiments on I^Cui-^Mg/Zn^G^ in 
the vicinity of the percolation threshold. 
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APPENDIX A: CALCULATION PROCEDURE 
FOR THE BOGOLIUBOV TRANSFORMATION 

We summarize below the steps required to obtain the 
eigenvector matrix for H, E satisfying the boson com- 
mutation relations (|2.39|) : 

• Diagonalizc H using the Lapack routines. This 
yields a set of eigenvalues u>i with corresponding 
cigen-subspaces generated by the eigenvectors 
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where n G {l,Pi} and where pi is the degeneracy 
of the eigenvalue and dimension of the subspace. 

• For the subspace i define E i: a 2N x p t matrix of 
the corresponding eigenvectors Ef . Form the block 
matrix M,- 



Mi = EjlEi, 

of size Pi x pi. 

• Invert M, to get M^ 1 . 

• Diagonalize M^~ , thus defining Ki and Df. 

Mr 1 = KiDiKr 1 . 

• Define the matrix: 



(Al) 



(A2) 



(A3) 



In this expression, the sign ± corresponds to the 
sign of uji . 

• Define the new matrix of eigenvectors for the sub- 
space i, Ei\ 



Ei — EiAi. 



(A4) 



• Repeat subspace by subspace to construct the 
eigenvector matrix E. 

APPENDIX B: STATISTICAL ERRORS 

this section we discuss the origin of the statistical er- 
rors we find from our numerical results. Consider, as an 
example, the lattice of size L 2 = 24 x 24. For x = 12%, 
we studied Mo = 620 different realizations of disorder. 
The average magnetization and root mean square (RMS) 
variation, AM, were found to be 



[M(x = 12%, L = 24)] P erc = 0.305 , AM = 0.0052 

(Bl) 

from which we estimate the error on the measure to be 
±%(i) — 



AM S 
VMo" 



(B2) 



In this example the estimated error is thus extremely 
small, around 0.1% and the errors rise to around 1% near 
the percolation threshold. This small error estimate is 
consistent with the statistical fluctuations observed in 
Figs. [Hand El 

For the example considered above, the ratio of the dis- 
persion to the mean value: 



AM 



[M] 



1.705%. 



(B3) 
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0.338 
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0.334 


0.0014 


0.419 


4 


0.329 


0.0022 


0.668 


6 


0.324 


0.0028 


0.864 


8 


0.319 


0.0036 


1.128 


10 


0.313 


0.0041 


1.310 


12 


0.305 


0.0052 


1.705 


14 


0.298 


0.0060 


2.013 


16 


0.290 


0.0070 


2.414 


18 


0.281 


0.0080 


2.847 


20 


0.272 


0.0096 


3.529 


22 


0.261 


0.0109 


4.176 


24 


0.249 


0.0135 


5.421 


26 


0.236 


0.0154 


6.525 


28 


0.221 


0.0183 


8.281 


30 


0.204 


0.0213 


10.440 


32 


0.186 


0.0237 


12.742 


34 


0.166 


0.0274 


16.506 


36 


0.147 


0.0291 


19.795 


38 


0.125 


0.0316 


25.280 


40 


0.109 


0.0330 


30.275 



TABLE I: Staggered magnetization for L = 24 x 24 and 
x £ [0,40]% 



L 


[M] pcrc 


AM 


10 


0.325 


0.0181 


12 


0.311 


0.0166 


14 


0.298 


0.0160 


16 


0.291 


0.0134 


18 


0.284 


0.0118 


20 


0.279 


0.0114 


22 


0.276 


0.0103 


24 


0.272 


0.0096 


26 


0.269 


0.0090 


28 


0.267 


0.0079 


30 


0.264 


0.0080 


32 


0.263 


0.0075 


34 


0.261 


0.0068 



pcrc 



TABLE II: Staggered magnetization for x = 20% and L € 
[10,34]. 



The ratio of the dispersion, AM, to mean value [M] perc , 
for fixed dilution L as a function of x and for fixed x 
and as a function of L are shown in Tables (JTJ) and (|iT|) . 
respectively. 

We can model this dispersion using three sources of 
variation: firstly, for a given x the number of magnetic 
sites varies from configuration to configuration. Sec- 
ondly, for fixed N the number of sites on the percolating 
cluster will also vary. Thirdly, there will also be a contri- 
bution from configurational fluctuations for a fixed num- 
ber of sites. We stress that all these contributions are 
quantum in origin. That is, the classical ground state is 
perfectly ordered for all concentrations above the perco- 
lation threshold, as discussed in the main text, hence at 
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the classical level, changing the number of sites, or local 
structures on the percolating cluster will not change the 
order parameter. However, the dilution reduces the local 
spin stiffness for spins in contact with non-magnetic sites 
and increases the zero point spin fluctuations. Hence 
these variations in number of sites and structure change 
the value of the order parameter. Indeed this point is 
already manifest by the fact that [M]p erc decreases with 
x. 

If iVj(L 2 , a;) is the number of sites for realization i, and 
the mean number of sites is defined in Eq. (|4.1[) then for 
the example considered we find 



N(L = 24,x = 12%) = 506.8 ± 7.7. 



Hence 



AN 



= 1.51%. 



(B4) 



(B5) 



To check the importance of fluctuations in the number 
of participating sites on the percolating cluster, (N perc )i 

we analyze the ratio (^p ctc ^ _ yje find 



N 



pcrc 



N 



= 0.999 ± 6.67 x 10~ s , 



(B6) 



which gives 



N 

N pcrc 



0.6% 



(B7) 



iV 



For a fixed number of magnetic sites, we can define the 
quantity 

AM 



as a measure of the configurational contribution to the 
dispersion in ground state order parameter values, where 
(• ■ ■ )„„„ is the disorder average over the restricted set 
of configurations with iVj = constant. For the example 
discussed here we find 



AM 



[M] 



pcrc 



= 0.4677% , 



(B8) 



from which we estimate the total dispersion 



/ AM 



\[M] 



P erc / tot 



AN 



x pcrc 

N 

■W] pcrc 



/ AM 

Xwr 



2% 



P erc / ma, 



N 



(B9) 

in good agreement with Eq. (|B|). The analysis can be 
generalized to the other values of x and L. For example 
for L 2 — 24 x 24 and x = 30% we have: 
AN 



N 
N nc 



N 



N, 



pcrc 



N 
AM 



[M] 



pcrc 



= 3.48% 



1.41% 



5.45% 



(B10) 



mag 



which correspond to 



AM 



10%, as obtained in 



perc / tot 

Table HJ Hence this analysis seems to account for the dis- 
persion in magnetization values to a good level of approx- 
imation. The three sources of dispersion are of the same 
order of magnitude as long as one remains well away from 
the percolation threshold. At low defect concentration 
(small x) it is the fluctuations in the number of magnetic 
sites that dominates. As x increases the fluctuations in- 
crease, as one might expect as one approaches the critical 
percolating regime, and at large x it is the configurational 
contribution for fixed particle number which dominates. 
At 40% dilution the dispersion in values approaches 30% 
of the mean order parameter value. Despite this large 
dispersion for this value of x the number of configura- 
tions, Ao = 2000, is large enough to keep the estimated 
error at the 1% level. 
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